function [vecTheta,Q,iter,exitflag] = Projection_LMoptimizer(vecTheta,setup,tolf,maxIter,lambdaInput,deltaOption,dispOn)
% This function minimizes the least squares objection function in projection
% using using the Levenberg Marquardt algorithm for the
% objective function min Q = sqrt(sum(f(x))^2)
jacobianVer     = 1;
gridOption      = 1; %1 for searching for Q_new<Q, 2 for optimizing Q_new
gridScaling     = [0.000001 0.001 0.1 0.3 0.5 0.8 1.0 1.2 1.5 3.0 5.0];
numVecTheta     = length(vecTheta);
f               = NLSobjectFunc(vecTheta,setup);
nf              = size(f,1);
Q               = f'*f;
iter            = 0;
exitflag        = 1;
updateJ         = 1;
maxNoDecrease   = 10;
countNoDecrease = 0;
diffQ           = 100;
epsValue        = 1D-6;
fp              = zeros(nf,numVecTheta);
fm              = zeros(nf,numVecTheta);
lambda          = lambdaInput;

delta           = ones(numVecTheta,1);

%parpool('local',setup.numCPUs); %start of multiprocessing
fastUpdateJ = 0;
while diffQ > tolf && iter < maxIter && countNoDecrease < maxNoDecrease
    iter = iter + 1;
    if updateJ == 1
        if fastUpdateJ == 1
            for i=1:numVecTheta
                if abs(delta(i))>epsValue
                    theta_p      = vecTheta;
                    theta_p(i,1) = vecTheta(i,1) + epsValue;
                    fp(:,i)      = NLSobjectFunc(theta_p,setup);
                else
                    fp(:,i)      = (J(:,i)*epsValue + f);   %ensures J(:,i) is unaffected
                end
            end
        else
            % The gradient
            for i=1:numVecTheta
                %parfor (i=1:numVecTheta,setup.numCPUs)   %for multiprocessing
                %parfor i=1:numVecTheta   %for multiprocessing
                % Positive step
                theta_p      = vecTheta;
                theta_p(i,1) = vecTheta(i,1) + epsValue;
                fp(:,i)      = NLSobjectFunc(theta_p,setup);
                if jacobianVer == 2
                    theta_m      = vecTheta;
                    theta_m(i,1) = vecTheta(i,1) - epsValue;
                    fm(:,i)      = NLSobjectFunc(theta_m,setup);
                end
            end
        end
        % Gradient
        if jacobianVer == 1
            J = (fp - repmat(f,1,numVecTheta))./epsValue;
        elseif jacobianVer == 2
            J = (fp - fm)/(2*epsValue);
        end
    end
    
    % Computing the new value of vecTheta
    if deltaOption == 1 %&& diffQ > 1D-4
        deltaOptionUsed = 1;
        JJ = J'*J;
        scalingCond = rcond(JJ + lambda*diag(diag(JJ)));
        if scalingCond > 1D-15
            delta = (JJ + lambda*diag(diag(JJ)))\J'*(-f);
        else
            delta = (J'*J + lambda*diag(diag(JJ)) + eye(numVecTheta)*1D-6)\J'*(-f);
        end
    else
        deltaOptionUsed = 0;
        delta = (J'*J + lambda*eye(numVecTheta))\J'*(-f);
    end
    
    
    % Evaluating fit at new value of vecTheta
    vecTheta_new  = vecTheta + delta;
    idxSearch     = 0;
    scaling       = 1;
    if any(isnan(vecTheta_new)) || any(isinf(vecTheta_new))
        exitflag     = 0;
        Q_new        = Q+1; % This ensures that we update lambda
    else
        if gridOption == 1
            f_new     = NLSobjectFunc(vecTheta_new,setup);
            Q_new     = f_new'*f_new;
            diffTempQ = Q_new - Q;
            if diffTempQ > 0 || isnan(diffTempQ)
                % We decrease the scaling to find a better point
                while (diffTempQ > 0 || isnan(diffTempQ)) && idxSearch < 10
                    idxSearch = idxSearch + 1;
                    scaling   = 0.8^idxSearch;
                    vecTheta_new = vecTheta +scaling*delta;
                    f_new     = NLSobjectFunc(vecTheta_new,setup);
                    Q_new     = f_new'*f_new;
                    diffTempQ = Q_new - Q;
                end
            else
                % We increase the scaling to boost the stepsize
                Q0 = Q_new;
                while diffTempQ < 0 && idxSearch < 10
                    idxSearch = idxSearch + 1;
                    scalingTmp   = 1.2^idxSearch;
                    vecTheta_newTmp = vecTheta +scalingTmp*delta;
                    f_newTmp     = NLSobjectFunc(vecTheta_newTmp,setup);
                    Q_newTmp     = f_newTmp'*f_newTmp;
                    diffTempQ    = Q_newTmp - Q0;
                    if diffTempQ < 0
                        scaling      = scalingTmp;
                        vecTheta_new = vecTheta_newTmp;
                        f_new        = f_newTmp;
                        Q_new        = Q_newTmp;
                        Q0           = Q_new;
                    end
                end
            end
        elseif gridOption == 2
            numGrid      = size(gridScaling,2);
            vecThetaGrid = zeros(numVecTheta,numGrid);
            fGrid        = zeros(nf,numGrid);
            QGrid        = zeros(1,numGrid);
            for i=1:numGrid
                vecThetaGrid(:,i) = vecTheta +gridScaling(1,i)*delta;
                fGrid(:,i)        = NLSobjectFunc(vecThetaGrid(:,i),setup);
                QGrid(1,i)        = fGrid(:,i)'*fGrid(:,i);
            end
            [~,idxSort]  = sort(QGrid);
            vecTheta_new = vecThetaGrid(:,idxSort(1));
            f_new        = fGrid(:,idxSort(1));
            Q_new        = QGrid(1,idxSort(1));
            scaling      = gridScaling(1,idxSort(1));
        end
    end
    DvecTheta             = abs(vecTheta_new  - vecTheta);
    [dThetaMax,theta1Pos] = max(DvecTheta);
    diffQ                 = abs(Q_new - Q);
    if dispOn == 1
        disp(['dQ          = ', num2str(Q_new - Q)]);
        disp(['Q           = ', num2str(Q_new)]);
        disp(['dThetaMax   = ', num2str(dThetaMax),' (at element ',num2str(theta1Pos),')']);
        disp(['lambda      = ', num2str(lambda)]);
        disp(['deltaOption = ', num2str(deltaOptionUsed)]);
        disp(['updateJ     = ', num2str(updateJ)]);
        disp(['Scaling     = ', num2str(scaling)]);
        disp(['iteration   = ', num2str(iter)]);
        disp(' ');
    end
    if scaling >= 1
        lambda = max(lambda/2,1D-20);
    else
        lambda = lambda*2;
    end
    if Q_new < Q
        vecTheta   = vecTheta_new;
        f       = f_new;
        Q       = Q_new;
        updateJ = 1;
        if fastUpdateJ == 0
            fastUpdateJ = fastUpdateJ + 1*0;%never use fast update of J
        elseif fastUpdateJ == 1
            fastUpdateJ = 0;
        end
        countNoDecrease = 0;
    else
        updateJ = 0;
        countNoDecrease = countNoDecrease + 1;
        lambda = lambda*5;
    end
    
    if dThetaMax < 1e-12
        if dispOn == 1
            disp('Change in x is less than 1e-12, so we stop the optimization');
        end
        diffQ = 0;
    end
end
if abs(diffQ) > tolf || iter > maxIter
    exitflag = 0;
end
%delete(gcp('nocreate')) %end of multiprocessing
if dispOn == 1
    disp('END OF OPTIMIZATION')
    disp(' ')
end
end
